
%Function to calculate consumer surplus of each project 
% and then compute average consumer surplus by country
function [cs , eV_dnk] = conSurplus(m, rho, par)

cs = zeros(1,4);

%First part stolen from model_constraints.m, convert rho to two matricies
% on for germany and one for denmark, each project by firm...
% Read root
rho_ger = rho(1:m.num_pro_ger * (m.num_firms_ger + 1) ,1);    
rho_ger_mat    = (reshape(rho_ger,m.num_firms_ger + 1,m.num_pro_ger))';

rho_dnk = rho(m.num_pro_ger * (m.num_firms_ger + 1)+1 : m.num_pro_ger * (m.num_firms_ger + 1) + m.num_pro_dnk * (m.num_firms_dnk + 1) ,1) ;    
rho_dnk_mat    = (reshape(rho_dnk,m.num_firms_dnk + 1,m.num_pro_dnk))';
    
%missing matrix part...

% 10.14.2013 - update for log state specification...
%extract the estimates....
fe      = par(1 : m.num_firms_ger);         % firm fixed effect (note: all the firms are active in Germany)
betaLogD = par(m.num_firms_ger + 1);              % coefficient on distance 
betaFor = par(m.num_firms_ger + 2);
betaState = par(m.num_firms_ger + 3);

%%
% Compute utilities: Germany

    %rho_ger_vec_f  = rho_ger_mat(:,1); %winning probabilities of fringe
    rho_ger_mat_nf = rho_ger_mat(:,2:m.num_firms_ger + 1); %win probs of active firms


    %contsturct mean portion of utilites for each firm and project 
    V_ger_nf  = -(repmat(fe',m.num_pro_ger,1) ...
                + betaLogD * m.dist_mat_ger ...
                + betaFor * m.foreign_mat_ger ...
                + betaState * m.state_mat_ger ...
                +  1 ./(1-rho_ger_mat_nf));  
    eV_ger = [ones(m.num_pro_ger,1) exp(V_ger_nf)];
    cs(1) = mean(log(sum(eV_ger,2)));
    cs(3) = (log(sum(eV_ger,2))'*m.ger_mw)/m.num_pro_ger; %Here we weight projects by megawatt capacity

% Compute utilities: Denmark

    %rho_ger_vec_f  = rho_ger_mat(:,1); %winning probabilities of fringe
    rho_dnk_mat_nf = rho_dnk_mat(:,2:m.num_firms_dnk + 1); %win probs of active firms


    %contsturct mean portion of utilites for each firm and project 
    V_dnk_nf  = -(repmat(fe(1:m.num_firms_dnk,1)',m.num_pro_dnk,1) ...
                + betaLogD * m.dist_mat_dnk ...
                + betaFor * m.foreign_mat_dnk ...
                + betaState * m.state_mat_dnk ...
                + 1 ./(1-rho_dnk_mat_nf));  
    eV_dnk = [ones(m.num_pro_dnk,1) exp(V_dnk_nf)];
    cs(2) = mean(log(sum(eV_dnk,2)));
    cs(4) = (log(sum(eV_dnk,2))'*m.dnk_mw)/m.num_pro_dnk; %Here we weight projects by megawatt capacity



